ANALYSIS OF THE TOOLKIT METHOD FOR THE 
TIME-DEPENDANT SCHRODINGER EQUATION 



Lucie Baudouin,Julicn Salomon, Gabriel Turinici 



Abstract. The goal of this paper is to provide an analysis of the "toolkit" 
method used in the numerical approximation of the timc-dcpcndent Schrodinger 
equation. The "toolkit" method is based on precomputation of elementary 
propagators and was seen to be very efficient in the optimal control frame- 
work. Our analysis shows that this method provides better results than the 
second order Strang operator splitting. In addition, we present two improve- 
ments of the method in the limit of low and large intensity control fields. 

1. Introduction 

The control of the evolution of molecular systems at the quantum level has been 
a long standing goal ever since the beginning of the laser technology. After an 
initial slowed down of the investigations in this area due to unsuccessful experi- 
ments, the realization that the problem can be recast and attacked with the tools 
of (optimal) control theory [TD] greatly contributed to the first positive experimen- 
tal results [TJ [Til H2 EE [20] . Ever since, the desire to understand theoretically how 
the laser acts to control the molecule lead the investigators to resort to numerical 
simulations which require repeated resolution of the Time Dependent Schrodinger 
Equation of the type (JTJ ; additional motivation comes from related contexts (online 
identification algorithms, learning algorithms, quantum computing [7], etc.). 

The numerical method used to solve the time dependent Schrodinger equation 
must provide accurate results without prohibitive computational cost. The conser- 
vation of the L 2 norm of the wave function ip(x,t) is also generally required for 
stability and as a mean of qualitative validation of the numerical solution. 

In this context, the second order Strang operator splitting is often considered [U 
[HI HHj- However, this method suffers from two drawbacks. First, the numerical 
error is proportional to the norm of the control which implies poor accuracy when 
dealing with large laser fields e(t) and make necessary the use of small time steps. 
Secondly, it requires at each time step three matrix products. This difficulty is 
enhanced in some particular settings e.g., in optimal control, where the matrices 
involved in the control term must be assembled online. 

Recently introduced, the "toolkit method" [2TJ [22] solves this last problem by 
precomputing a set of elementary matrices, used in the numerical resolution. Each 
matrix is associated to (one or several) field values and enables to solve the evolution 
over one time step. This algorithm has been used in various frameworks and shows 
excellent results. It has also been coupled successfully with optimal control and 
identification issues [3]. The dependence on the L°°-norm of the control, which is 
a restriction of the Strang method, is also improved by the "toolkit method" as it 
will be shown in our analysis. 

The goal of the paper is to provide a (first) numerical analysis of the "toolkit 
method" . Our mathematical tools are related to that in [4] (but for a different 

l 



ANALYSIS OF THE TOOLKIT METHOD FOR THE TIME-DEPENDANT SCHRODINGER EQUATION 



setting; see also [9] [16] for connected results); the treatment here is different because 
of the quantization appearing in the values of the control e(t) which impacts both 
the mathematical analysis and the numerical efficiency of the method. The analysis 
enables us to propose two possible improvements. 

The paper is organized as follows: after having introduced the model and some 
notations in Section [21 the toolkit method is presented and analyzed in Section [3] 
An improvement of this method in the limit of small control fields is introduced 
in Section 2] A second improvement, in the limit of large control fields is given in 
Section [5] Finally, Section [5] gathers some numerical results. 



In this section, we present the Schrodinger Equation that will be considered in 
the paper and some useful notations. 

We consider the time dependent Schrodinger equation (TDSE) 



This equation governs the evolution of a quantum system, described by its wave 
function ip, that interacts with a laser pulse of amplitude e, the control variable. 
The factor /i is the dipole moment operator of the system. The Hamiltonian of 
the system is Hq = —A x + V where A x is the Laplacian operator over the space 
variables and V — V(x) the electrostatic potential in which the system evolves. 
We refer to [13] for more details about models involved in quantum control. Note 
that to obtain Eq. ([T]), one has considered the laser effect as a perturbative term, 
so that the control term s(t)fi(x) is obtained through a first order approximation 
with respect to s(t). While often considered, this approximation fails at describing 
some models involving non linear laser-dipole interaction, see e.g. [15] . Conse- 
quently, the norm of the field cannot be always considered as a small parameter, 
and numerical solvers have to tolerate large controls, as the one described here after. 

Throughout this paper, T > is the time of control of a quantum system. The 
space L p (0,T; X), with p 6 [l,+oo) denotes the usual Lebesgue space taking its 
values in a Banach space X. The notation 14 Al,1 (0,T) corresponds to the space of 
time dependent functions belonging to L 1 (0,T;R) such that their first time deriv- 
ative also belongs to L^O, T). We denote by L 2 the space L 2 (M 3 , C) and by VF 2 '°° 
and H 2 the Sobolev spaces Fy 2 -°°(R 3 , R) and # 2 (R 3 ,C). The space £{H 2 ) is the 
space of linear functionals on H 2 . One can refer to [5] (or the introduction of [6]) 
for more details about the definitions of these functional spaces. 

Finally, in order to introduce some numerical solver of |T]), let us consider an 
integer N and At > such that NAt = T . We introduce the time discretization 
(tj)o<j<N of [0, T] with tj — jAt and we also denote by tj + i the intermediate time 



2. Model and notations 



(1) 




U + f )At. 



2 



Let us first recall some basic results of existence and regularity of the solution of 
the TDSE. These are corollaries of a general result on time dependent Hamiltonians 
(see [M], p. 285, Theorem X.70). 
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Lemma 1. Let fi G C{H 2 ), V € W 2 ^ , e G £ 2 (0,T) and ip G H 2 . The 
Schrodinger equation 

M f i^(i) = (ffo-^(*))#) M 3 x(0,T) 

lj W(0)=^, K 3 , 

ftas a urape solution tp € L°° (0, T]H 2 )n W 1,oo (0, T; L 2 ) suc/i i/wrf 

||^(*)IU°°(0,T;iP) + ||9t^(*)|U«( 0) r;ia) < C|H| £ (jj2) ||e|| W 1 ' 1 (0,T) I! tpO \\ ■ 

Moreover, for all t G [0,T], = ||Vo||^- 

It is also well known (see [6] for instance) that for any T > and do £ H 2 , if we 
have e(t) = e G R, independent of time t, the Schrodinger equation 

f id t <f>(t) = (H - t ie)cj>(t) ) K 3 x(0,T) 
1 0(0) = K 3 

has a unique solution = S(i)0 o such that G C([0,T];iJ 2 ) n C 1 ([0, T];L 2 ), 
where (5(t))teK denotes the one-parameter semi-group generated by the operator 
H a - fie. Moreover, for all t G [0,T], S(t) G £(# 2 ) and we have 

S(t)0 o eC(0,T;tf 2 ), V0 o e^ 2 ; 

II^WIU(^) < l + CT<if, VtG[Q,r];A' = A'(|H| jC (fl a) ,e„ u «); 

(3) ||5(t)0 Q || i2 = \\Ml*, V0 o eL 2 ,Vie]R; 
S(0) = Id, 

S(t + s)=S(t)S(s), Vs,<eR. 

Therefore, the solution of Eq. |[2j) is obtained equivalently as a solution to the 
integral equation 

ip(t) = £(t)t/>o + i / 5(i - s)(e(s) - e)fJtip(s) ds. 
Jo 

3. The toolkit method 

We now present the toolkit method and describe the corresponding error analysis. 

3.1. Algorithm. In this method, we assume that the control field e satisfies the 
following hypothesis: 

(H) Vt G [0,T], e(t) G [ 

The values of the control field are discretized according to: 

(4) et = e min + £Ae, I = 0...m, 

with m = gmax A ~ £min . Here, the values eg have been here uniformly chosen in the 
interval [e m i n , £max] ■ If some properties of the field are known, e.g. its mean value 
or its variance, some improvement of the method can be obtained by optimizing 
the distribution of the values eg. More generally, this topic enters the field of scalar 
quantization, that will not be considered in this paper. We refer to [17j and the 
references therein for a review of standard methods in this domain. In order to 
solve numerically equation @, the toolkit method proceeds as follows. 
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Algorithm 1. (toolkit method) 

(1) Preprocessing. Precompute the "toolkit", i.e. the set of propagators: 

S t (At) for I = (),■■■ ,m, 

where (Sg(t)) t £K. denotes the one-parameter semi-group generated by the 
operator Hq — fisi, the sequence (e^)^=o,- - ,m, being defined by ((4]). 

(2) Given a control field e 6 L 2 satisfying El and %/jq — ip$, the sequence 
(V , i f r )j=o,...,./v that approximates (4>(tj))j=o,...,N, is obtained recursively by 
iterating the following loop: 

(a) Find: 

lj = argmin^ =1) ... )m {|£(t i+ i) -e e \}, 

(b) Set^f +l = S tj {At)^f . 



In this "toolkit" approximation, we consider that the changes in the Hamiltonian 
H(t) := H Q — fis(t) can be neglected over a time step At. In this way, if Ae = 
(infinite toolkit), and for a relevant time discretization, the simulation correspond- 
ing to piecewise constant control fields is exact. Such a property does not hold 
with methods that approximate the exponential, e.g. the second order Strang op- 
erator splitting. Indeed, these approaches introduce an algebraic error, due to the 
non-commutation of the operators Hq and [i that is consequently proportional to 
||e]|i«>(o,T)- 

Remark 1. In the original form of the toolkit method pi 1211 [2"2"] . the mid-point 
choice proposed in Stev \2a\ of Algorithm [71 is not considered. Yet, the introduction 
of this strategy enables us to improve the order of the method (see the analysis 
hereafter). 

3.2. Analysis of the method. Let us now present an error analysis of the toolkit 
method. More precisely, this section aims at proving the following result: 



Theorem 2. Let e € Vl /2 ' oo (0, T) and ip the corresponding solution of @. Let ip K 
be the approximation of ip obtained with Algorithm^ Given At > and Ae > 0, 
there exists Ai > 0, A2 > 0, that do not depend on ||£||l°°(o,t) such that: 

(5) U(T)-ij K (T)\\ L 2 <X 1 Ae + X 2 At 2 . 

Moreover, there exists v\ > 0, v 2 > depending on ||£||h' 1 ' 1 (o.t) such that: 

(6) \\ip(T) — ^ K {T)\\ H 2 < v\Ae + v 2 At 2 . 

Remark 2. This result shows that the toolkit method enables to work with large 
control fields, transferring the computational effort due to such cases to the pre- 
processing step: given Ae, the computational cost of this step only depends on the 
norm ||£||l~(o,t); i- e -, on Hypothesis Wu 



Proof. To obtain ([5]) and ©, we will focus first on the local error, i.e. the approx- 
imation obtained on one time step [tj, tj+i). 

The sequence )j=o,...,iV is a time discretization of the solution of 

f id t ib K (t) = (H - f ie(t))^ K (t), K 3 x(0,T) 
( ' \ ip K (0) = il> , R 3 
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where the space variable has been omitted and s(t) — eg. is constant over each inter- 
val fe,t i+ i[= \jAt,(j + l)At[, withj = 0,...,iV-l. We denote by (Sj(t)) j=Qj N _ x 
(instead of Sg , ) the one-parameter semi-group generated by the operator Ho — fj,S£. 
and we introduce S(t) = e(t) — e where e (instead of eg ) is the constant value of 
e{t) over [tj,tj+i\. Therefore, the solution ip of $2$ is actually the solution of the 
integral equation, settled for t 6 [tj,< J+ i[: 

(8) ip{t) = Sj(t - tj)ip(tj) + i [ Sj(t- s)n5(s)ip(s) ds. 

Jtj 

For the upcoming calculations, one should notice that we have 

(9) W j+ i)\ < ^ 
and that for all t € 

(10) \S(t)\ <i(A£ + ||e|| i . {0 , T) At). 
We consider the following decomposition: 

i>{T) - ip K (T) = </>(T) - Sy_i(Af)V(tw_i) 

N-2 
3=0 

+ S N -i(At)...S (At)i> -'<p K (T) 
where the last line is equal to since ip K satisfies © on [0,T]. 

From now on and in all the following sections, we will consider either that 
II^oIIl 2 = 1 or that HV'ollff 2 = 1- From (|3|), we know that the operators Sj are 
isometrics in L 2 . Therefore, the use of a triangular inequality brings 

JV-l 

(11) u(t) - 4> k {t)\\ l , < £ sjiAtwmiv . 

3=0 

We will thus calculate and estimate in L 2 -norm for all j the difference 

Tp(tj + i) — Sj(At)ip(tj) — i Sj(tj + i — s)5(s)/j,ip(s) ds 

tj+i 

Sj(tj + i — s)<5(s)/i (ip(s) — Sj(s — tj)ip(tj)) ds 

i 

rtj+i 

i / Sj(tj + i — s)S(s)fiSj(s — tj)4>(tj) ds. 

Sj(tj + i — s)S(s)/j, (ip(s) — Sj(s — tj)ip(tj)) ds 

i 

[■tj + i 

(12) + i I S(s)(p j (s)ip(t j )ds 

Jtj 

where tpj(s) := Sj{tj + \ — s)jiSj(s — tj) G C{L 2 ). 
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In what follows, we work in parallel on L 2 and 7J 2 -estimates of ip(tj+i) — 
Sj(At)ijj(tj). We will need basic L 2 and i? 2 -estimates of ip{t) — Sj(t — tj)ij){tj) 
for the study of the first integral term of (fT!?|) . the second one will be dealt with 
using a Taylor expansion of 8(t). 

From Lemma[TJ (jHJ and ([lOII , it is easy to obtain coarse estimates of ip(t) — Sj (t — 
tj)ip(tj). Indeed, for all t in [tj,tj + {\, one can write 



L 2 



< 



i / Sj(t — s)/x(5(s)-0(s) ds 



\Sj(t — s)fj,S(s)i/j(s)\\ L2 ds 



< AtM^- (Ae + \\i\\ L ^ T) At) Uo\\ L 2 



< 



\n\\c(V) (Ae + p|| r oo (0 ,T)At) At 



(13) 

and the if 2 -estimate gives 

(14) U(t)-S j (t-t j )il>(t j )\\ Ha < K\\e\\ w i,i (0tT )Mc(m) (Ae + ||e|U- (0iT) At) At 

where K = K{\\ii\\ c ^ H 2^ e max ) is a generic constant that estimates every H^H/j/jm. 
Therefore, we can obtain more accurate estimates of the first integral term of (p~2|) . 
Thanks to (fT3|) . we obtain 



L 2 



i / Sj(tj + i — s)<5(s)/x (tp(s) — Sj(s — tj)ip(tj)) ds 



1 



< dHU(£») (Ae+||e|U» (0iT) At) At sup \\4>(s) - Sj(s - t^fe)IU 2 
A te[t 3 ,t J+1 ] 

< i|MlWA £ +|| £ |U~ (0 ,T)At) 2 At 2 
(15) < ^llMllV) (Ae 2 + p||ioo ( o, r )Ai a ) At 2 . 

Working now on the ff 2 -estimate, we deduce from ([14")) in the same way that 



(16) 



Sj(tj+i - s)6(s)fj,(ip(s) - Sj(s - tj)ij){tj)) ds 



L- 



< ^K\\plU{ H -)\\4w^{q,t) (Ae 2 + ||e||^ (0iT) At 2 j At 2 . 



In the two cases (L 2 and H 2 ), estimates are stronger than the ones we look for, 
and we can focus on the second integral term of (fT2f we want to deal with. 



We first consider 

(17) vr> Mi+i] -^£( H2 ) 

s i-> Sj{t j+1 - s)fiSj(s -tj) 
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and note that for all ijj e H 2 , \\(pj(s)t/j\\H2 = ||Sj(ij+i — s)fiSj(s — tj)ip\\ H 2 < 
K\\ip\\ff2 so that 



(18) 



Vs e [tj,t j+1 ], \\<Pj{s)\\c(H») < K\\n\\ C (spy 



Let us now consider the derivatives of ipj{s). Since (Sj(t)) te s. denotes the one- 
parameter semi-group generated by the operator Hq — jii, the £(H 2 ) identity 

d t S {t) = -iiHo-^Sjit) 

holds and minor calculations give, Vs 6 [tj,tj+\], 

9 s (Pj(s) = iSj(tj +1 - s)[H ,fj]Sj(s -tj) 
dssfjis) = SjOj+i- s)[[H ,fi},H - iie]Sj(s-tj). 

Therefore, 

(19) |]cWiOOIU(i^) < K\\[H ,»]\\ C(H2) 
\\ds S <pj(s)\\£( H z) < K\\[[Ho,h],Hq -^e]\\ c{H2) . 

If we consider the L 2 -analysis of the method, then (pj(s) € C(L 2 ) and Vs e [tj, tj+i], 

\\<Pj(s)\\c{L*) < IIm]U(£ 2 ) 

(20) \\d.M*)\\cW < \\[H , f A\\c(L^ 
\\9s S <fj{s)\\c(L^ < \\[[Ho,l4,Ho- (J>e]\\ c(L2 y 

Let us now write the third order Taylor expansion of t i — * 5(t) = e(t) — e in a 
neighborhood of tj + i: 

6(a) = S(t j+ + (s-t j+ i)S(t j+i ) + ^(s-t j+i ) 2 S\e(s)) 
= 6(t j+ i) + (s - t J+i )i(t j+ x) + ha - t J+h ) 2 e{B{s)), 



with 8(s) 6 [tj , tj+i] . We now focus on estimating the term i 
By means of (f^Uj) and the L 2 -norm conservation, we obtain 



5(s)(pj(s)ip(tj) ds. 



<*(VfiVj(s)V'(*j) ds 



L 2 



< - \\n\\c(L*) AeAt, 



s - t j+ i )e[t 



V 1 



- f'j+i) / s + s ) - 





0)^(t,-)df 



(*'+*) 



< Y^I|[^0,/-i]||z:(L2)||£||L°=(tj,t 3 + 1 ) 



sd u Lp(u)ip{tj) duds 
At 3 



L' 2 



ANALYSIS OF THE TOOLKIT METHOD FOR THE TIME-DEPENDANT SCHRODINGER EQUATION 



and 



ftj+i 1 ( v 2 



< ^7llMll£(L^)||£||L~(t 3 ,t 3 + 1 ) Ai3 - 



24 

L 2 



Combining these results with (|15p. we estimate (fT2")l as follows: 
IWi+J-SjiAtWm* < IMc(lz) (A e 2 Ai 2 + ||e||^Ai 4 ) 

+ 1 (2||[if ,M]]]£(^)P||oc + ||/i||£(^)||£']|oo) At 3 , 

with || • ||l=(o.t) = II ' lloo- By means of (fTTj) . the global L 2 -estimate is then: 

JV-l 

U(T) - ^ K (T)\\ L2 < Yl llV^i+i) - ^-(At)^fe)|| L2 

< (Ae 2 Ai + pHLAt 3 ) + | ||M||^)A£ 

+ J (2||[^o,M]IU(^)PIU + IHk(£»)l|e||<») A ^ 

and |5]| can be deduced with the following constants Ai and A2 independent of 

I|£||l~(o,t) : 

Al - ||Ml£(L2)(l + A£Ai|H| £(L2) ), 

A2 = - IMI^PII^At + — ||[Ho,M]||£(z»)l|e||oo + |jlMl£(L2)||e||oo. 

Let us now prove the H 2 estimate. By means of (|T5|) . (|T5|) and p9[) and keeping 
in mind that K is a generic constant depending on ||/^||£(h 2 ) and e max , we can 
repeat the previous analysis to find the local estimate: 

H{tj+i) - SjiAtWtj)]]^ < K\\ii\\jnH*)\\e\\ w i.no,T) (Ae 2 Ai 2 + ||e||i«- (0 ,r) A * 4 ) 

+ KAeAt + K(\\lH , f i}\\ c{H 2 } \\e\\ L ~ + \\i\\ L ») At 3 . 



Since one can prove that we can actually write a more precise estimate of Sj(At) 
and replace K by 1 + CAt (see properties ©), we get: 

\\S 3 (At)\\ c{H 2 } < 1 + CAt. 

and since we have the following intermediate result, where M > depends on 
IImI|£(£T 2 )! £ max an d T but is independent of N: 

JV-l 

Y J (^ + CAt) N - j At < M, 
3=0 
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The global estimate is obtained as follows: 

AT-l 

U{T)-^ K {T)\\ H , < ^K^-^Uitj+^-S^AtMtj)^ 

3=0 
AT-l 

< 2 (1 + CAt) Jv ^|| e || W r l .i (0l D (Ae 2 Ai 2 + Hell^At 4 

JV— 1 

+ ]T (1 + CAt)»-* (As At + (|| [H , h]\\c(h 2 ) ¥\ 

3=0 

< M||e|| WI . 1( o,T) (Ae 2 Ai + Hell^At 3 ) + MAs 
+ M(||[flo,M]IU(fl»)l|e||oo + P1U) At 2 . 

We finally get v\ and v% and conclude the proof of Theorem [5] 
j/i = M(l + ||e|| w i,i (0ir) rAeAt) 

^ 2 = M (||e|| M ri.i(o l r)||e||LAt+||[flo,A*]||/:(fl»)l|e||<» + l|e||oo). 



□ 



Remark 3. XTie estimate ([5]) zs consistent with the fact that Algorithmic used with 
a relevant time discretization is exact for the piecewise constant control fields. 



4. Improvement in the limit of low intensities 

We now describe a way to improve the time order of the previous algorithm. 
Since some constants in the following analysis depend in this case of the L°°-norm 
of the field and the method requires that the toolkit size scales A£ 3 (e max — e m i n ), 
it applies in the case of (L°°-) small control fields. 

4.1. Algorithm. The algorithm we propose mixes the toolkit and the splitting ap- 
proaches, in the sense that it applies sequentially various operators to correct the 
third order local error that appears in the proof of Theorem [2] 



Algorithm 2. (Improved toolkit method for low intensities) 

(1) Preprocessing. Precompute the "toolkit", i.e. the set of propagators: 

S e (At) for£ = 0, ••• ,m, 

where (Se(t))teR denotes the one-parameter group generated by the operator 
Hq — iaei, the sequence (se)e—o t ... , m , being defined by (f4]). Include in this 
set the two special elements: 

and the initial exponents ag and (3$ such that (e being extended as an even 
function on [-T,0]): 

e(At) - e(0) ., . . ?x 
a := M =^\) + 0{At% 

e(ti)-2e(ii)+e(0) 
00 := A^ -=e(t h ) + 0(At 2 ). 
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(2) Given a control field e € L°° satisfying^ and ip 1 ^ = fl a °Q^°ijjo, the 
sequence (ipj K )j=o,...,N that approximates (ip(tj))j=o,...,N, is obtained re- 
cursively by iterating the following loop: 

(a) Find: 

tj = argmin f=1) ... !m {|e(^ + i/ 2 ) - e e \}, 

(b) Compute ctj and (3j such that: 

(21) a, := £ ^- £ fe) =e{t J+k ) + 0{M% 

e(t j+1 )-2e(t i+ i)+e(tj) 

(22) j := -AlL! K -^=e(t lH ) + 0(Ati). 



(c) Set = Se^AtJO^e^^ 



IK 
3 ' 



In many cases, e.g. in the experimental frameworks, only the values of the field 
can be handled. The use of exact values for the time derivatives has then to be 
avoided when possible. This motivates the introduction of approximations (|21[) and 
([22]) of s(tj) and s(tj) in the latest definitions. The analysis presented hereafter 
shows that this does not deteriorate the order of the method. 

In this method, one must perform two online matrices exponentiations. By 
working in a basis where one of these two matrices is diagonal, the cost of Step dc| 
can be reduced to one exponentiation, making the cost of this method equivalent 
the second order Strang operator splitting. 



4.2. Analysis of the method. We can now repeat the analysis that has been 
done in the proof of Theorem [2] to obtain the following estimate. 



Theorem 3. Let e G W 2,oo (0, T), ip be the corresponding solution of @ andip IK 
the approximation of tp obtained with Algorithm^ Given At > and Ae > 0, 
there exists X[ > 0, X' 2 > 0, with X[ independent of ||£||l°°(o,t) such that: 

mT)~^ IK (T)\\ L2 <X[Ae + X' 2 At 3 . 

Proof. In the framework of this new algorithm, we note that on every time interval 
}tj,tj + i[, the approximation ip IK is the solution of the evolution equation: 



(23) 



id t ip IK (t) = {H - ne)V K (t), K 3 x (tj,t j+1 ) 

Tp IK (t+) =n a ie^i/j IK (tj) r 3 



where we set ip(0 ) = ipo- We will keep the notations (Sj, S(t), <p,...) of the proof of 
Theorem [21 and we first focus on the local error analysis. We consider the following 
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decomposition: 



tf>(T) - ^ IK (T) = if,(T) - S N -i(M)Sl a > r - 1 Q fl ' r - 1 il>(t t r-i) 



N-2 



+ s N - 1 (At)n aN - i e' 3N ~ i . . . s j+1 (At)n a ^+ i e^+ i 

3=0 

x (^(*i+i) - Sf(At)n°'eftv(*j)) 
+ s N ^ 1 {At)n aN - i Q l3N - 1 . . . s o (At)n ao e 0o ipo - ^ ik (t) 

where the last line is equal to since ip IK satisfies ^ on [0,T]. 

The operators Sj are isometries in L 2 , we will consider that HV'oIIl 2 = 1 and we 

also have, for all j 

(24) n^e^ = e 3[*b,/*]A* 8 e &/*At a = w+ fe[# 0jM ] + l -^i?J At 3 + ld 0(At 6 ) 



and thus 
(25) 

Therefore, the use of a triangular inequality brings 



n^e^lL™ < i + o{At 3 ). 



N-l 



(26) U(T) - ^ IK {T)\\ L , < (1 + <D(At 2 )) ^ |^fe + i) - ^-(At^e^foOL 

and we will calculate and estimate in L 2 -norm for all j the difference 
ip(t j+1 ) - SjiAt^e^ipitj) 

= SjiAt^itj) +i f 3+1 s 3 {t J+1 - s )s(s)^(s) ds - s,-(At)n a *e"*v(*i) 
it, 

= Sj(At) (Id - 0^6*) V(ij) + * I Sj{t 3+1 - s)5(s)(j,ip{s) ds. 

Jtj 

We define Y(s) = ip(s) - S 3 (s - tj)Q a J 8& ip(tj) for all s e [tj,tj+i] and obtain 

(27) Y(t j+1 ) = Sj(At) (Id - fi^G^) i)(tj) 

+ i ( 1+ S J (t J+1 -s)S(s)piY(s)ds + i ( 1+ S(s)ipj(s)n^Q^ip(tj)ds 



where fj(s) := Sj(tj+± — s)fiSj(s — tj) and its derivatives have been estimated in 
L 2 in (12"01) . As we did in Theorem [5J we start with an estimate of the first integral 
term of l|27p. For all t £]tj,tj+i], we can write: 

Y{t) = i>(t) - Sj(t - t^itj) + Sj(t - t )^{t 3 ) - Sj(t - tj^e^ipitj) 

Sj{t - s)S(s)fiip{s) ds + Sj{t - tj) (Id - il^Q^) V(ij)- 



Moreover, for all t &}tj, we have 



\\Y(t)h* < 



Sj(t — s)S(s)^iip( s ) ds 



+ \\Sj(t- tj) (id - n°*e&) tp(tj)\ 



L- 
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The operators Sj are isometries in L 2 and Vi e [0, T], HV'WIU 2 = HV'oIIl 2 - There- 
fore, we deduce from (|24|) that 



s 3 (t - 1 3 ) (id - n«ie^) ^fe)|| L2 < ( ^ ||[ffo,M]IU( i2) + IWU(l>) ) Ai 3 +o(A; 6 ) 



24 



Since it is clear that we also have 



S^(t — s)S(s)fiip(s) ds 



< - 



(Ae+||e||i-, ( o, r) At) At, 



L 2 



one can finally deduce that: 



*3 + l 

/ / Sj(tj+i — s)6(s)fj,Y(s) ds 



L 2 



1 



< -||/i||£(L 2 )(Ae+p|| L ^ (0 ,T)At)Ai sup ||F(t)|| L 2 



(28) 



< 



tG[tj,tj_(_i] 

^IHl£ ( L 2 ) (Ae + ||e|| ioc( o,T)At) 2 At 2 + 0(AeAi 4 ) + 0(At 5 ). 



We focus now on the first and third terms of (|27p . Using (|24ll . we get 



•S«A/)(I«l- <-r'B'\) ,-(/,) = -.S(A/i ( ig[ff , M ] + ^) r(/,)Ar-:- ( .(/ y iC7|A/"). 



Let us then consider the second integral term of (|27|) . On the one hand, we consider 
the fourth order expansion of 6 — e — e in a neighborhood of 

= «5(t i+ i) + (s t j+i )6(t j+i ) + \{s t jH ) 2 S(t j+h ) + i(, - tj . + i )V 3 )(#( S )) 
= <J(t i+ i) + (s - * i+ i)e(t i+ i) + i(« - t j+j _) 2 e(t j+i ) + l(s- t J+h )V 3 > (9(s)) 



2' 6 ' 



where 0(s) S [tj, On the other hand, we calculate and/or estimate the four 

corresponding terms in 



tj+i 



i / <5(s)p i (*)n a *e"<v(*i)d»- 
it, 

From © and (|2"0j) . the term of order gives: 



i / % + i)^( a )n°*e'V(* J > 



< 



,AsAt. 



L 2 
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For the term of order 1, we can write 



i (s- t i+ i)£(^. + i)^-( s )o^e^vfe) ^ 

i e {t J+i ) ^ s(w(t j+ i + s) - Vj (t j+ i - s))n a i&^{t ) ds 
i : ( / ; . i j / / " ^ sd^jiu^Q^ipitj) duds 



i s(t 



■3+1 



ft . . 1 —s 



s(d u ipj(tj) + (u- t j )T(u))ft a ie^'i/3(tj) duds 



£{tj ^S 3 (At) [Fo,M]^e*V(*i)At 3 

■ iAt r t , 1+s 



+ 1 



Jo Jt i-s 

J +2 



s(u - tj)r(u - ^)fi Qj e ft ip(tj) duds 



= ZLSj (At) [H , ^(t 3 )At 3 + Sj (At) [H , u]^p(t j )0(At t 



12 



r At rt J+i +. 



' t . . 1 — s 



i(u - ^)t(m - tj-)fi^e^^(tj) duds 



where we used dH}, (glj and (gD]) and the function r : s G [0, At] t (s) G C(L 2 ) 
is defined as the function that appears in the following expansion of d u (fj around 
tj, for any ip G L 2 

d u (pj(u)ip = d u tpj(tj)i/)+(u-tj)T(u-tj)<if) 

= iSj(At)[H , fj]ip + (u- tj)T(u - tj)ip. 

Using the estimate (coming from (|20|)) 



(29) 



\\ t ( s )\\c(l 2 ) < [H ,n],H - as 



C(L 2 ) 



Vs G [0,At], 



along with HVK^OIU 2 = HV'oIIl 2 = 1, (EH an d J2S|) we find that for all j, 



< (a 3 +0(At 2 )) 



iAt ,* J+ i+ S 



4 At rt, +i +« 



i(it - tj)T(u - t j )0, a ^e^tlj(t j ) duds 



L 2 



JO Jt.i-s 
3 + ? 



s(u-tj) [Ho,n],H - ue \\ip(tj)\\ L2 duds. 



< 



24 



[H ,/j],H - fie 



C(L 2 ) 



At 4 + 0(At 5 



We also prove easily that for all j, 

\\S 3 (At) [Ho,u]^j(t 3 )0(At e )\\ L2 = O (At 6 ) . 

For the term of order 2, using (f22l) . ([24ll and (|20p and the first order expansion 
of ipj around tj, (fj(s)ip = ip 3 (tj)ip + (s — tj)9(s — tj)ip for all tp G L 2 , defining 
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9 : s e [0, At] i ► 0(s) e C(L 2 ), we can write 

l J* 3+1 l -{ S -t 3+h fs{t ]+h )^ 3 {s)^Q^^t 3 ) ds 
= ie(t, + i) / -(s - t i+h ) 2 (tpjitj) + (s - tj)9(s - tj)) n^Q^^tj) ds 



ie(t + i) 

3 2 -Sj (Ai)^ aj © ft V'(^)At 3 



24 

ie(t i+ i) f'i+i 



(s - t j+ i ) 2 (s - t.j)o{ s - tj)n a 'e^ip(tj) ds 



(At) ^(t 3 )At 3 + S, (At) ^(t 3 )0(At 6 ) 



ie(t j+ i) r*j+i 



(s - t J+h f(s - tj)B{a - tj^QtoWi) ds. 



Using (|2"0|) , we get the estimate ||0(s)||,c(£2s < H ,fi] 
using it with ([22]) and ([25]) . we obtain that for all j, 



Vs 6 [0,Ai], and 



L 2 



< -(/^ + 0(Ai 2 )) / (.s-t J+ i) 2 fe- S )||[i/o,M]|| £(L2) ll^)ll L2 ^ 



< 



< 



-(P J+ O(At 2 ))\\[H ^}\\ c{L2) 



ft 



At 






;t-) 







|8 ll^o./i]!!^ At 4 + 0(At 5 ). 

and we also prove easily that ||5,- (At) ^i/ , (ii)C(At 6 )|| L2 = O (At 6 ) . 
Combining these results with (|28|) into equation (|27| . we obtain: 



< 



Sj(At) (id - n°*e^) vfe) + i / 5(s)^-(s)o^e /3 ^(^) ds 



L 2 



Sj(tj + i — s)S(s)/jY(s) ds 



< 



2 ^ £ At+|| 



L 2 

[#0, /•*],#() - 



£(L 2 ) 



At 4 



+ § H[ffo,dlL (La ) At 4 + ill^ll 2 :^) (Ae 2 At 2 + \\i 
+ 0(AeAt 4 ) + 0(At 5 ) 



2 A/ 4 
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We have now a local in time estimate that should be traduced in a global one, and 



from (126 



we get 



r K (T)\\ 



< 



< 



{l + 0(At 2 )) 



N-l 



3=0 



J 1 

48 



Ae 



ip(t j+ i) -5j(At)n^e*v(* 

otjT 



\L 1 



+ 0(AeAt 3 ) 
The result follows, with 



24 

\\[H ,^}\\ c{L2) At 
0(At 4 ). 



[H ,/j],H 
T 



lie 



At 3 



+ 2^11^) (A^At + PIII^^At 3 ) 



-I 

2 1 



Mlli (i 2)(l+A £ Ai) 



and 



24 



£(L=) 



48 



II [#o, 



T 



3 (0,T) 
□ 



In this theorem, the constants A' 2 depends on ||e||_L°°(o,T) through the commutator 



[H ,fi\,H - fie 



that appears in 



This contrasts with the result obtained in 



Theorem [3] The explanation of this situation comes from the fact that the norms 
of ifij(s) := Sj(tj + i — s)/j,Sj(s — tj) (defined in (|T7|> ) and its first derivative does 
not depend on ||£||l=°(o,t), whereas its second derivative does. Thus, errors in 
Algorithm [5] depend on L°°-norm of the control field as in the case of the second 
order Strang operator splitting. Although these two methods present the same 
computational complexity, the order of Algorithm [5] is higher when Ae scales At 3 . 

5. Improvement in the limit of large intensities 

We now describe a way to improve the time order of the Algorithm [T] in the case 
of large intensities. The following method enables to replace Ae by AeAt in the 
estimates. 

5.1. Algorithm. The algorithm we propose improve the accuracy in the approxi- 
mation of e. This improvement is obtained by using two toolkit elements instead 
of one at each time step. 



Algorithm 3. (Improved toolkit method for large intensities) 

(1) Preprocessing. F 'recompute the "toolkit", i.e. the set of propagators: 

Si(At) fori = (),••• ,m, 

where (Si(t))t<£R denotes the one-parameter group generated by the operator 
Hq — \iEi) the sequence (ei)i=o.... . m , being defined by (j3|). 

(2) Given a control field e G L°° satisfying 1771 and ipQ K — ipo, the sequence 

{4>j K )j=o,--.,N that approximates (ijj(tj))j = o n, is obtained recursively by 

iterating the following loop: 
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(a) Find £j such that: 

(b) Compute ctj and (3j such that: 

otjE^ + (3jEt ]+ i = e(t j+1/2 ) 
(30) aj+Pj = 1 



(c) Set tj+i = S ij+1 [Atfi S tj {M) a ^ 



JK 
j ' 



In this method, one must perform two online matrices exponentiations. The 
cost of the corresponding step, namely Step He] can be reduced to three matrix 
products when precomputing the mappings between the diagonalization basis of 
two consecutive toolkit elements. 

Remark 4. Another way to reduce the cost of this step, is to quantify the val- 
ues of otj (and fjj) and precompute a toolkit containing elements of the form : 
Si j+ i(At)Pj (At) a i . This method is tested in Sec. [3] 

5.2. Analysis of the method. We can now repeat the analysis that has been 
done in the proof of Theorem [2] to obtain the following estimate. 

Theorem 4. Let e € W 3 '°°(0, T), ip be the corresponding solution of ^ and ip JK 
the approximation of tp obtained with Algorithm [J] Given At > and Ae > 0. 
there exists X" > 0, X' 2 ' > 0, both independent of \\£\\l°°(o.t) such that: 

||V(T) - t JK (T)\\ L 2 < X'lAsAt + X'^At 2 . 

Proof. In this algorithm, two control fields are involved successively in the propa- 
gation over the interval [tj,tj+\]. As in the previous proofs, we introduce 8(s) = 
e(s) — e(s), with 

si. s e [tj,tj + ttj-Atf, 

e £j+1 s S [tj + a 3 At,t j+1 [. 
Note first that for all s £ [tj, tj+i] 

1 

At 



s(s) 



(31) |<y( a )|<Ae + -||e|U-( ,r) 

and denote by (Sj(t))^_ N _ 1 and (S'j(t)j^_ Q N _ 1 the one-parameter semi- 
groups generated by the operators H Q — ii£g j and H — /if^.+i respectively. 
Following the same analysis as for Algorithm [TJ we set {ipj )j=o,„,,N as t ne time 
discretization of the solution of: 

' id t ^ JK {t) = (H - fie(t)) ^ JK {t), R 3 x (0, T) 
</r>*(0) = Vo, M 3 



(32) 



where e(t) (defined right above) is constant over each interval [tj,tj + (XjAt[ and 
[tj + ctjAt, tj+\[, with j = 0, N — 1. In the same way as we obtained ([5]), the 
solution if) of ^ satisfies, 

ptj+otj At 

ip(tj + ajAt) — Sj(ajAt)ip(tj) —if Sj(tj + ajAt — s)fj,5(s)ip(s) ds 
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and 

ip(t j+ i) = S'jiPjAtfipitj + otjAt) - i / S'j(t j+1 - s)ti5(s)ijj(s) ds. 

Jtj+ctjAt 



As in (jl2|) , it gives rise to: 

S'jitj+i - s)/j,5(s)ip(s) ds 

J tj-\-a.j At 
ptj-\-ajAt 

Sj{PjAt)Sj(tj + ctjAt - s)fiS(s)i!(s) ds 
S'^tj+t - s)fiS(s) (ip(s) - S'jis - tj - ctjA^ipitj + ajAt)) ds 

Itj+ajAt 

/tj+OLj At 
S'0jAt)Sj{tj + ctjAt - s)[i5{s) (ip(s) - Sj(s - tj)if)(tj)) ds 

+i / 5(s)ip' J (s)S J (a j At)i}(t j ) ds 

Jtj+ajAt 

ftj-\-OLjAt 

(33) +i S'jiPjAQSWvjisMtj) ds. 

where cp'j (s) := S'j(tj+i—s)fj,Sj(s—tj — ajAt) and^-(s) := Sj(tj +a>jAt~s)nSj(s — 
tj). As in the proof of Theorem [2] (see right above (jTTJ)) we use the appropriate 
decomposition 

V>(T) - ^ JK (T) = ip(T) - S'^iPN^A^SN-iiaN^At^itN-x) 

N-2 

+ S' N _ 1 (p N -iAt)S N -t(a N -- i At) . . . S' j+1 (p j+ iAt)Sj+x(a j+ iAt) 

x(^%+i) - S'j{PjAt)Sj{ajAt)^{tj)) 
+ S , N _ 1 (PN-iAt)S N -i(a lf -iAt) . . . S' (p Q At)So(a At)^ - yj JK (T) 

where the last line is equal to since ip JK satisfies (J32j) on [0,T]. We have the 
corresponding estimate (see (fTTj) ) 

iV-l 

U(T)-t JK (T)\\ L 2 < W^tj+J-S'jiPjAVSjiajAtMtj)]]^ 

3=0 

and we will thus calculate and estimate in L 2 -norm for all j the four terms of (|33| . 
As in (JT5| . but using now the new estimate (|3Tj) of 5, the two first terms of the 
right hand side of can be respectively estimated by: 



(34) 



S'j(t j+ i - s)fi6(s) (ip(s) - Sj(s - tj - ctjAt)i){tj + ctjAt)) ds 

tj+OLj At 



L- 

1 I A -2 a 4.2 , 1 \\±\\2 A ,4\ 



< PiMc(L') (A^A^ + ^m\u {Q>T) At 
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and 



(35) 



S'^A^Sjitj + aj At - s)^5(s) - S^s - t 3 )^{t 3 )) 



L' 1 



1 



2x, ' J ' 5plli- (0 ,T)A* 4 



Let us now focus on the third and fourth terms of (|33|) . We have: 



5(a)ip' j (a)S j (a j At)ip(t j ) ds = 



tj + At 



6(s)<p'Jtj + a j At)S j (a j At)ip(t j ) ds 



tj -\-OtjAt 



d u $j( u ) du Sj(ajAt)i()(tj) ds 



5{s) 

tj-\-OLj At J tj-\~OtjAt 

S(s)S' j (P J At)^S j {a j At)ilj(t j ) ds 

tj-\-OLj At 



5( S ) 



tj+ajAt Jtj+ajAt 



d u ip'Au) du Sj(ajAt)ip(tj) ds 



and 

ptj +aj At 



S'(^-Af)<5(s)^ i (s)V'(^) ds 



tj-\-a.j At 



5(s)S' j (j3 j At)<p j (tj + ajAt)il}{t j ) ds 



tj -\~Otj At 



S^S'^At) 



tj-\-ctj At 



d u ipj(u) du ip(tj) ds 



ti+ctjAt 



5(s)S!j (fy At)nSj (ajAt)^ ) ds 



tj -\~OLj At 



5(s)S'MAt) 



tj-\-ajAt 



d u ipj(u) du ipitj) ds. 



By means of (f3"0)l . we have: 



6(s) ds 



e(s) - s(t j+1/2 ) ds+ e(t j+1/2 ) - e(s) ds 



e(s) - e(t J+ i/ 2 ) rfs 



where 9(s) S [ij, ij+i]. Consequently, 



(36) 



5(s)S'J(3 j At)iJ,S j (a j At)ip(t j ) ds 



L- 



< — ||Mll£(^)l|e||L~(o,T)Ai 3 
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From (EH and (J3TJ) , we obtain 

rtj+i 



(37) 



5(a) 



tj+ctjAt 



tj+ajAt 



duf'Au) du Sj(ajAt)ip(tj) ds 



L- 



< \(%\\[H ,lA\\c{L*) fAe+i||e|U- (0ir) At) At 2 



and similarly, we find that 

ftj-i-atj At 



(38) 



tj+otjAt 

S{s)S' J (p j At) I d u <pj{u) du ipfe) ds 



L- 



< \<%W*AUw (Ae+i||e|| £ c O(0 ,T)At) At 2 



Combining (J34J) , (J35J) , (36]), (J3TJ) and (38]), we obtain: 



\^(t j+1 ) - S'APjAVSjiajAt)^)] 



L- 



< 



\l(L2) 



(As 2 + i||e||i O o (0 , r) A* 2 ) At 2 



+ 



24' 



|£(L 2 )I|£||L°°(0,T) 



At 3 



+^||[ffo,//]|U(i=) fAe+i||e|| £ oo ( o lT )A^ At 2 



The global estimate follows 



JV-l 



U(T) - ^ JK (T)\\ L2 < ll^fe+i) - S'^jAt^iajAt)^ 



3=0 



< 



1 



A^ + -||£ 



2 \\ a \\L°°(0,T) 1 



At 2 TAt 



+ ^Mc(L*)\\e\\L~( ,T)TAt 2 

+ ^\\[Ho,fl]\\c(L2) (Ae+ i||£|| £ ~ (0 ,T)^ 



At TAt 



and the proof of Theorem [4] is complete, with 

1, 



A" = i,\\[Ho^]\\c(L^T + M 2 c(L2) TAe, 



A 2 = \\\[H ,[4\\c(l2)\\£\\l°°(o,t)T+ ^M\c(l--)\\s\\l~(^t)T 

+ \M\l { L*)¥\\L~{Q,T)TAt. 



□ 



6. Numerical results 

In this section, we check numerically that the order of the estimates we have 
obtained in this paper are optimal, and we compare computational costs of the 
methods. 
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6.1. Model. In order to test the performance of the algorithms on a realistic case, 
a model already treated in the literature has been considered. The system is a mol- 
ecule of HCN modeled as a rigid rotator. We refer the reader to [T5] for numerical 
details concerning this system. 

As a control field, we use an arbitrary field of the form e(t) = e max sin(wi), with 
E max = 5.10 -5 and u> — 5.10~ 6 . The parameters are chosen in accordance with 
usual scales considered for this model. The use of an analytic formula for the field 
enables us to work with exact values, i.e. to test the cases Ae = 0. 

6.2. Orders of convergence. To test the time order, we first work with Ae — 0, 
with various values of At. The numerical orders correspond to the ones obtained 
in our analysis. Curves of convergence are depicted in Fig. [TJ 



-e — Splitting a=1.9979 




FIGURE 1. Error with respect to At, when Ae = for toolkit 
method and Improved toolkit I method, and whenAe = cAi for 
Improved toolkit II method. Here, ip num stands for the approx- 
imation of ip when using the toolkit method, the second order 
Strang operator splitting, the Improved toolkit I method and the 
Improved toolkit II. The coefficient a is the regression coefficient. 

The order with respect to Ae is also obtained numerically by using a small 
time step. In this test, the numerical order is consistent with the one obtained in 
Theorem [5J The convergence with respect to this parameter is presented in Fig. [5] 

6.3. Computational cost. In a second test, we compare the computational costs 
of the methods. To do this, we look for the values of N = and m — that 
enable to reach a fixed arbitrary error of Tol — 5.10 -3 (recall that in any case the 
error cannot exceed 2). For sake of simplicity, we only test powers of 2. In this test, 
we also include the quantified version of the Improved toolkit II which is described 
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I 



-Toolkit a=0. 98935 



10~* 
10' 



10 

Ae/e max 



Figure 2. Error with respect to Ae, when At is small. Here, 
tjj num stands for the approximation of ip when using the toolkit 
method. 

in Remark [5] In our case, the parameters a and /3 were quantified among 100 values 
uniformly distributed in [0, 1]. 







Matrix products 


777 = 


Strang Op. Splitting 


16384 


32768 




Toolkit 


8192 


8192 


16384 


Improved toolkit I 


1024 


2048 


16384 


Improved toolkit II 


4096 


12288 


16 


Quantified Improved toolkit II 


4096 


4096 


6400 



Table 1. Values of numerical parameters corresponding to a tol- 
erance error of Tol = 5.10 -3 . 



These tests show that toolkit methods always give better results as the second 
order Strang operator splitting. 

The two improvements we propose in this paper enable to reduce respectively the 
global number of matrix products and the size of the toolkit, which is in agree- 
ment with the analysis we have done. Note that the second improvement reduce 
significantly preprocessing step. This fact makes feasible the quantified version of 
it, which requires intrinsically a larger toolkit. 
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